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Abstract 

Background: When studying metabolism at the organ level, a major challenge is to understand the matter 
exchanges between the input and output components of the system. For example, in nutrition, biochemical models 
have been developed to study the metabolism of the mammary gland in relation to the synthesis of milk components. 
These models were designed to account for the quantitative constraints observed on inputs and outputs of the 
system. In these models, a compatible flux distribution is first selected. Alternatively, an infinite family of compatible set 
of flux rates may have to be studied when the constraints raised by observations are insufficient to identify a single flux 
distribution. The precursors of output nutrients are traced back with analyses similar to the computation of yield rates. 
However, the computation of the quantitative contributions of precursors may lack precision, mainly because some 
precursors are involved in the composition of several nutrients and because some metabolites are cycled in loops. 

Results: We formally modeled the quantitative allocation of input nutrients among the branches of the metabolic 
network (AIO). It corresponds to yield information which, if standardized across all the outputs of the system, allows a 
precise quantitative understanding of their precursors. By solving nonlinear optimization problems, we introduced a 
method to study the variability of AIO coefficients when parsing the space of flux distributions that are compatible 
with both model stoichiometry and experimental data. Applied to a model of the metabolism of the mammary gland, 
our method made it possible to distinguish the effects of different nutritional treatments, although it cannot be 
proved that the mammary gland optimizes a specific linear combination of flux variables, including those based on 
energy. Altogether, our study indicated that the mammary gland possesses considerable metabolic flexibility. 

Conclusion: Our method enables to study the variability of a metabolic network with respect to efficiency (i.e. yield 
rates). It allows a quantitative comparison of the respective contributions of precursors to the production of a set of 
nutrients by a metabolic network, regardless of the choice of the flux distribution within the different branches of the 
network. 
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Background 

When studying metabolism, it is important to elucidate 
how fluxes are distributed among the different pathways 
of the metabolic network with respect to the available 
quantitative information about the system behavior. Sev- 
eral methods can be used to address this issue. The 
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first approach consists of building a mechanistic descrip- 
tion of transformations and identifying the regulations 
involved in the system. Continuous dynamical models are 
often used for this purpose, especially when time-series 
responses to different treatments are available to infer the 
dynamics of the network. Static approaches such as Petri 
net can also identify qualitative distributions of fluxes in 
a metabolic network [1], and their stochastic extensions 
can even take into account stoichiometric and kinetic 
information [2]. With a complementary approach, one 
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can study a system at the functional level, based on the 
study of fluxes at steady states. This is the purpose of 
the Flux Balance Analysis (FBA) framework, which has 
evolved considerably over the past decades [3,4]. With 
an FBA, the identification of external regulations is not 
necessary because it is assumed that the global behav- 
ior of the system can be modeled by optimizing linear 
combinations of selected fluxes (i.e. the objective func- 
tion). Roughly, the methods developed in this field aim 
to explore a convex space of plausible flux distributions 
and to study the extreme flux distributions obtained when 
optimizing a linear objective function. It allows check- 
ing whether an extreme flux distribution is consistent 
with experimental data and to predict new experimen- 
tal observations [5-8]. Nonetheless, the consistency of 
the solutions obtained by FBA depends on the quality 
of the constraints integrated in the model. To overcome 
this limitation, several extensions have been proposed 
in the literature. These extensions can be broken down 
into two parts. The first incorporates additional biolog- 
ical knowledge, such as reaction thermodynamics [9] or 
multioptimization [10]. The second is based on the use of 
FBA to globally analyze large-scale metabolic networks. 
For example, in flux variability analysis [11], the minimum 
and the maximum flux for each reaction in the network 
are computed under some (sub-)optimal conditions. 

In this paper, our main purpose is to extend this frame- 
work to study the variability of a metabolic network at 
the level of efficiencies instead of fluxes. Indeed, when 
studying a metabolism at the organ level, a major chal- 
lenge consists in comparing the efficiency, or yield rates, 
of two metabolic situations, i.e. the response to vari- 
ous input patterns. A typical example of such studies are 
those concerning animal nutrition, which aim to predict 
the quality and quantity of animal production (meat, fat, 
milk, etc.) in response to breeding factors. In this field, 
the energy and protein conversion efficiencies are derived 
from the study of the flux distribution of input nutrients 
between the different branches of the metabolic network 
[12]. More generally, although its definition depends on 
the field of application, the concept of efficiency is often 
linked to energy, mass growth and protein conversion 
[13-15]. However, the computing of efficiencies is prone 
to difficulties and errors, since it requires computing the 
quantities precursors sets which are required to explain 
the measured composition of outputs, although some of 
the internal products may be recycled within cycles [16]. 
With this goal, we defined the allocation of an input 
towards an output (with the abbreviation AIO) to be the 
proportion of a matter component (such as carbon) in a 
given input flux that is recovered in the selected output 
flux. Our definition is based on the choice of a material 
component, such as carbon or nitrogen, which allows a 
comparison of the contributions of input metabolites to 



the composition of output products. It can be seen as 
a yield rate, which is uniform among all the outputs of 
the system, allowing a precise understanding of the pre- 
cursors used. As a first methodological contribution, we 
prove that AIO can be uniquely described and computed, 
even in the case when there are metabolic system cycles, 
by a matrix whose coefficients are nonlinear functions of 
the flux variables. Introducing nonlinearity in the defini- 
tion of AIO cannot be avoided because of the presence of 
these cycles. 

Studying the variability of AIO within a complete space 
of plausible flux distributions requires the solving of non- 
linear optimization problems which are underdetermined 
in tangible applications. As a second methodological con- 
tribution, we have proposed efficient algorithms to com- 
pute lower and upper bounds for AIOs over the fam- 
ily of flux distributions which are compatible with both 
the system's stoichiometry and the experimental datasets, 
regardless of the choice of a flux distribution for the 
internal branches of the network. An important aspect is 
that when the metabolic network is provided with input- 
output data, the complete space of plausible distributions 
appears to have a relatively small size, and can therefore 
be studied with our method. Our framework is depicted 
in Figure 1. 

Our main example of application is related to milk pro- 
duction. In this context, several models have been intro- 
duced, in relation with the aforementioned classification 
of models. One class of small-size dynamical mechanis- 
tic models predicts the blood flow and input nutrients 
of the metabolic system (i.e. the mammary gland) [17,18] 
or, alternatively, the nutrients produced by the metabolic 
system (in terms of milk composition) [19,20]. Another 
class of models predicts the distribution (i.e. partitioning) 
of fluxes over the pathways from the input and output 
nutrients of the system [21,22]. In the latter family, both 
dynamic and static approaches exist. Indeed, numeri- 
cal models based on mass-action equations were initially 
proposed to describe the fluxes to and from individual 
metabolites, with different levels of description [21,22]. 
To that goal, optimizers were used to determine a reason- 
able set of parameter estimates for the dynamical model of 
the system. Although such a set of parameter estimates is 
not unique. With a totally different approach, in another 
study, the main biochemical pathways in the different 
metabolisms were integrated in a generic stoichiometric 
model called the metabolic spreadsheet [15,23]. Its struc- 
ture was based on a restricted number of intermediary 
metabolites called carbon-chain pivots, which correspond 
to important cross-over points between metabolic path- 
ways. This allowed the computation of the flux rates for all 
reactions, constrained by a general rule stating that there 
is no accumulation of intermediary metabolite of cofac- 
tors. A study of manual calculations performed with the 
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Figure 1 Main functionalities of the analysis workflow. First, extremal vertices of the simplex polyhedron of plausible flux distributions have to 
be computed, including the case when this space is not bounded. Then, formal algebra is required to obtain a symbolic representation of AIO 
matrices, expressed as a formal function, where the variables are the coefficients of a plausible flux distribution. Finally, extrema of the AIO 
coefficients are computed among the complete simplex of plausible flux distributions, either with the existing optimization routine or with 
dedicated local-search algorithms. 



metabolic spreadsheet showed that it works, in practice, 
by maximizing an objective function (ATP production) in 
a convex solution space [24]. Therefore, this model can be 
considered as an application of the Flux Balance Analysis 
framework [25]. 

Nonetheless, in this field of study, there has been little 
discussion on the impact of the choice of a single model 
among several (possibly infinitely many) reasonable mod- 
els [14]. To investigate this issue in a more automatic 
way, we compared the aforementioned conventional mod- 
els to the convex space of plausible flux distributions 
associated with steady states of a model of mammary 
gland metabolism, as computed in FBA. We checked 
the consistency of extreme flux distributions of nutrients 
with experimental data of the mammary gland and milk 
production, including the contribution of nutrient input- 
output and isotope balance studies [26,27]. Based on our 
AIO computation framework, our analysis highlighted 
that the metabolic behavior of the mammary gland cannot 
be modeled by maximizing ATP production or by opti- 
mizing a linear combination of flux variables of the model 
of the mammary gland. In other words, although an infi- 
nite number of flux distributions are compatible with the 
data, none are extreme within the space of feasible mod- 
els. Selecting any of these nonoptimal flux distributions 
to predict the system behavior appears difficult without 
additional experimentation. 



To gain better understanding of the system response 
regardless of the choice of a flux distribution for the inter- 
nal branches of the network, we applied our method to 
estimate the variability of AIO coefficients in our model 
and compared the effects of two different diets on mam- 
mary gland metabolism. Our results suggest that the 
bounds of AIO are sufficient to distinguish the effects of 
different nutritional treatments without selecting a flux 
distribution for the internal reactions of the metabolic 
network by any method - optimization of a linear combi- 
nation of fluxes or a residual score. Overall, the complete 
study suggests considerable flexibility in mammary gland 
metabolism. It provides a view of the functioning of the 
system although its internal processes still cannot be clar- 
ified because of limitations on experimentation on large 
animals such as ruminants. 

Results 

We first investigated the set of flux distributions that are 
compatible with the stoichiometry of our mammary gland 
model (depicted in Figure 2), without taking datasets into 
account. The model exhibited a large variability, since 
thousands of extreme pathways could be identified. 

We then successively computed the set of flux distribu- 
tions compatible with the model and the real datasets of 
lactation metabolism in dairy cows. The datasets are given 
in Table 1. They include a control diet (Ctrl), a diet related 
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Figure 2 Simplified view for the stoichiometric model of ruminant mammary metabolism with no long-chain fatty acid oxidation. The 

complete model is detailed in a SBML Additional file 1 . Eleven nodes are considered as pivots (green nodes), that is, intermediary metabolites which 
are not accumulated in the cell. The model is built by adding 26 specific reactions to the ruminant mammary to the generic model of [23]. A 
treatment is characterized by a set of input nodes (yellow nodes), quantitatively described in mmol/h/half udder. Nutrients contained in the milk are 
considered as output nodes. The node "Fatty acid synthesis" is an abstraction for 14 reactions corresponding to C(4:0), C(6:0), ...C(16:0) from 
primerC2 and primerC4. Note that, depending on the treatment, the role of nonessential amino acids may change from input to output according to 
the balance of the amino-acid considered. The node "peptide" summarizes the ATP cost of protein synthesis and protein degradation. Notably, the 
stoichiometry of reactions was adjusted in order to balance carbon exchanges, including C02.This is a key point in order to compute exact allocation 
tables in the following stages. Additional literature-based information allows us to generate additional linear constraints on some reactions fluxes. 



to an increased protein supply through casein infusion 
into the duodenum (CN), and a complementary dataset 
(HB) previously used in a mechanistic model of the mam- 
mary metabolism [21]. For the three datasets, the set of 
plausible flux distributions - solutions to a linear system 
detailed in Eq. (1) - was an unbounded convex cone of 
dimension 5. In all cases, the solution spaces shared the 
same set of five independent variables. This suggests that 
there are five independent levels of variability within the 
system: peptide hydrolysis (7? 64 ), NADPH oxidation (#19), 
OAA^PYR (R u ), OAA^G3P {R 15 ) and G3P^G6P (R s ). 
Therefore, by including this dataset the model became 
a fairly small and constrained network. Nevertheless, it 
was not uniquely determined since there were still several 
degrees of freedom. 

Investigating the relevance of the optimization strategies 
for mammary metabolism 

The balance between the ATP generated by the system 
and the ATP used by the system (including the ATP cost 
of milk component synthesis) was computed for the three 
datasets (Ctrl), (CN) and (HB). The results are detailed in 
Table 2. In this table, the FBA approach based on opti- 
mization of the ATP balance is compared to the natural 
functioning of the system. 



First, we considered the manual computation of fluxes 
with a tool named "metabolic spreadsheet" [15,23]. We 
applied the rules of no accumulation on any intermedi- 
ary metabolites and cofactors with a utilization rejecting 
all the cycles. As expected, since some cycles are ATP- 
consuming, both flux distributions were equivalent in 
both approaches (Table 2). The ATP balances were respec- 
tively 3081, 2045 and 6628 mmol/h/half udder (i.e. 318 
mol/d/udder) in the (Ctrl), (CN) and (HB) datasets. 

According to this model, the ATP generated is estimated 
to be 6500 mmol/h/half udder (i.e. 312 mol/d/udder) 
while 2125 mmol/h/half udder (i.e. 102 mol/d/udder) 
were estimated to be used for milk component synthesis 
[31]. Therefore, the ATP-balance obtained for this model 
was slightly smaller than the optimal one obtained for our 
mammary gland model. 

Independently of the approach considered and the 
model or dataset at hand, the remaining ATP (ATP bal- 
ance), after use for milk synthesis, appeared to be rather 
high and nonxconstant. Indeed, as shown in [31], an ATP 
balance of 1250 mmol/h/half udder (60 mol/d/udder) can 
be expected to be used for other functions not accounted 
for in the models, such as maintaining membrane poten- 
tial and synthesizing nucleic acids. The numerical val- 
ues obtained here were far from this estimated balance, 
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Table 1 Net uptake and milk component output of the mammary gland in three treatments 

Input or output flux (Ctrl) [28,29] (CN) [28,29] (HB)[21] 

mmol/h/ mmol/h/ mol/d/ mmol/h/ 







half udder 


half udder 


udder 


half udder 






Glucose input "> 


237 


232 


12.21 


254 




Vg 5 


Glycerol input 


5.84 


5.74 


0.033 


0.69 




V96 


Acetate input 


510 


462 


18.42 


384 




Vg 7 


BHBA input (2) 


84 


167 


7.25 


151 




vgs 


Lactate input 


0 


0 


0.023 


0.48 




V62 


3C(n:m)-acycoA+glycerol- 
3P^ triglyceride'- 3 '' output <4) 


32.96 


39.11 


1.52 


31.67 




Fatty acid output (synthezized)' 5 ' 












VlOO 


C(4:0) 


10.08 


11.59 


0.46 


9.48 




V101 


C(6:0) 


4.51 


5.58 


0.18 


3.79 




V\ 02 


C(8:0) 


2.23 


2.87 


0.10 


2.06 




Vl03 


C(10:0) 


4.66 


6.46 


0.19 


3.96 




Vl04 


C(12:0) 


4.23 


6.12 


0.17 


3.56 




Vl05 


C(14:0) 


13.90 


17.89 


0.45 


9.31 




Vl06 


C(16:0) 


18.82 


21.44 


0.64 


13.40 




Vgg 


Lactose output 


73.80 


83.52 


3.81 


79.28 




Amino acids balance' 61 i.e. entry or output 












Vl 28 


Alanine input 


3.11 


0 


0.105 


2.19 


Alanine catabolism 


V121 


Alanine output 


0 


3.26 


0 


0 


Alanine synthesis 


Vl 19 


Arginine input 


4.40 


4.48 


0.526 


10.96 


Arginine catabolism 


Vl34 


Asparagine output 


0 


0 


0.023 


0.48 


Asparagine synthesis 


V125 


Aspartate output 


3.43 


4.13 


0.247 


5.15 


Aspartate synthesis 


V 122 


Glutamate output 


0.54 


6.33 


0.230 


4.79 


Glutamate synthesis 


V131 


Clutamine input 


1.22 


1.79 


0.072 


1.50 


Glutamine catabolism 


V]20 


Glycine output 


4.98 


3.44 


0.248 


5.17 


Glycine synthesis 


VM4 


Proline output 


10.65 


10.99 


0.670 


13.96 


Proline synthesis 


V136 


Serine output <7) 


7.21 


7.50 


0.090 


1.88 


Serine synthesis - Serine 
used in other pathways 


V118 


Histidine input 


0.23 


0 


0 


0 


Histidine catabolism 


V113 


Isoleucine input 


2.19 


3.57 


1.518 


31.63 


Isoleucine catabolism 


VU4 


Leucine input 


2.02 


3.76 


0 


0 


Leucine catabolism 


vim 


Lysine input 


2.68 


3.58 


0.191 


3.98 


Lysine catabolism 


Vl 1 1 


Threonine input 


0.35 


0 


0 


0 


Threonine catabolism 


Vl 1 5 


Valine input 


2.54 


3.86 


0.438 


9.13 


Valine catabolism 


V107 


Peptide output (8) 


124.5 


150.0 


7.2 


149.17 




Additional constraints 












V82 


NADPH through ICDH pathways (9) 


30% 


30% 


30% 


30% 




,/ 0-7,, 
V58 = 03V82 


NADPH through Pentose Phosphate (9) 


70% 


70% 


70% 


70% 




V56 = 3V 62 


C(n:m)->C(n:m)-acylCoA 


98.87 


1 1 7.32 


4.56 


95.00 




FA primer from Acetate' 10 ' 












v 53 


C(4:0) 


0 


0 


0 


0 




V54 = Vgo 


C(6:0) 


2.256 


2.790 


0.091 


1.896 




V55 = Vgy 


C(8:0) 


1.113 


1.437 


0.050 


1.031 
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Table 1 Net uptake and milk component output of the mammary gland in three treatments (Continued) 



V51 = V S2 


C(10:0) 


2.331 


3.230 


0.095 


1.979 


VS6 = Vg 2 


C(12:0) 


2.116 


3.061 


0.086 


1.781 


V S7 = Vg 3 


C(1 4:0) 


6.951 


8.946 


0.223 


4.655 


V88 = V94 


C(1 6:0) 


9.410 


10.719 


0.322 


6.698 


Other constraints' 111 










V 24 


Lactate ->■ Pyruvate 


0 


0 






V44 


Alanine catabolism 




0 






V76 


Alanine synthesis 


0 




0 


0 


V83 


Asparagine synthesis 


0 


0 






V40 


Histidine catabolism 




0 


0 


0 


V36 


Leucine catabolism 






0 


0 


V33 


Threonine catabolism 




0 


0 


0 



Data are renormalized in mmol/h/half udder, (a) Control diet (Ctrl) [28,29] (b) Higher protein supply by casein infusion in the duodenum (CN) [28,29] (c) Generic 
dataset (HB) [21]. This table is used to parameterize input and output vectors v ( and v 0 together with additional biological linear constraints on some reaction fluxes. 
1 1nput i.e. taken up by the stoichiometric system considered (i.e.net uptake in our example for the mammary gland). 
2 /?-Hydroxybutyrate. 

3 Total triglycerides secreted in milk considering that milk fat was composed of 1 00% triglycerides and that all the triglycerides were secreted in milk fat [21 ]. 
4 Output i.e. leaving the system (secreted in milk). 

5 All fatty acids synthesized within the mammary gland i.e. all C4 to CI 4 and 50% of CI 6 [21]. 

6 The balance between amino acid net uptake and amino acid net output in milk protein is calculated with established rules [21,30]. If the balance is positive it 
corresponds to an amino acid input that will be catabolized. If this balance is negative, the values corresponded to an amino acid output (that will be synthesized to 
meet its requirement for milk protein). 

7 Serine output corresponded to Serine synthesized minus Serine utilized in other pathways i.e. Serine required in addition to Ser uptake to synthesize milk protein. 
8 Peptide output: number of peptide links required to synthesize the proteins exported out of the system (i.e. in milk protein). 
'Hanigan, 1994 [21]. 

10 Fatty acid primers were synthesized for 50% from acetate and for 50% from BHBA except C4 FA primer which was supposed to be synthetized only from BHBA [21]. 
11 Set at zero because their inputs or outputs are set at zero (to avoid futile cycle). 



Table 2 Three different computations of ATP balance for 
the mammary gland in different treatments 


Biological model 
Criteria of selection of a 
solution 


Dataset 


ATP balance 


Proteic 
turnover 


Mammary-gland model 
(Figure 3) 


(HB) 


6628 


0 


Manual study based on 
cycle removal [23] 


(Ctrl) 


3081 


0 




(CN) 


2045 


0 


Mammary-gland model 
(Figure 3) 


(HB) 


6628 


0 


ATP optimization 


(Ctrl) 


3081 


0 




(CN) 


2045 


0 


Model of Hannigan- 
Baldwin [21] 


(HB) 


4375 = 6500-2125 


0 


Study of numerical 


(Ctrl) 


Non available 





equilibria of an ODE 
model [31] 

(CN) Non available 

Three natural assumptions are considered to model the mammary gland 
behavior: removal of all cycles, optimization of ATP production and study of the 
equilibria of a dedicated ODE-based model. All models exhibit considerable 
variability in their ATP balance (in mmol/h/half udder), which contradicts the 
assumption about the behavior of this organ. Moreover, quantitatively, the 
computed ATP balances are much higher than recent measurements. This 
suggests that ATP maximization cannot be considered as a natural objective 
function to model cow mammary behavior. 



suggesting that each of these three models have a nonrel- 
evant ATP balance. 

In addition, both in the ODE model and the ATP- 
optimization approach, peptide hydrolysis was obtained 
at zero, implying an absence of any protein turnover. This 
contradicts all the observations about this pathway: con- 
siderable use of this pathway has been evidenced in several 
publications, although the peptide hydrolysis rates dif- 
fered significantly depending on the technique used for 
the measurements: peptide hydrolysis (7? 64 i.e. mammary 
protein degradation) spans from 0.25, 0.23 to 0.67 of pep- 
tide synthesis (7?63 i-e. total mammary protein synthesis) 
in (Ctrl), (CN) and (HB), respectively [29,31]. 

Overall, we concluded with this analysis that the energy- 
based optimization function may not allow an appro- 
priate simulation of mammary gland metabolism. This 
was expected, considering that the system is studied at 
the complete organ level, involving competing processes 
which can rarely be modeled with a single linear objective 
function. 

Exploring all extreme flux distributions in a refined simplex 

In order to study the variability within the space of plau- 
sible flux distributions and to identify alternative rele- 
vant optimization strategies, an additional constraint was 
placed on the ATP balance of the system. We considered 
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that an ATP-balance of 1250 mmol/h/half udder (60 
mol/d/udder) was a relevant measure for this study [31], 
although we checked that all the results discussed below 
were still valid when introducing a 10% tolerance on 
this estimation of ATP-balance. Providing a bound for 
ATP yielded a new constraint on several variables of 
the system. For the two datasets (Ctrl) and (CN), the 
flux G3P^G6P (Rs) was no longer an independent vari- 
able in the system. The simplex, which was previously 
unbounded, appeared to be bounded with four inde- 
pendent variables: OAA^PYR(7?i 4 ), OAA^G3P (R 15 ), 
NADPH oxidation (7? 19 ), peptide hydrolysis (7? 64 ). 

Applying flux variability tools [32], it appeared that 
the minimum of each of the fluxes V14, V15, V19, V(,i was 
equal to zero for all treatments. The maxima of the fluxes 
(V14, V15, V19, v 64 ) were (1831, 1831, 669, 305) for the (Crtl) 
treatment and (795,795,22,133) for the (CN) treatment 
This suggests that the (Ctrl) treatment generates a more 
flexible space of plausible flux distributions than the (CN) 
treatment. 

We computed all extreme vertices of the simplex of 
plausible flux distributions for the two treatments (Ctrl), 
(CN). The simplex structure of the polyhedron implies 
that the optimum of any linear combination of metabolic 
fluxes involved in the model is either uniquely attained 
for one of these extreme points or attained by all flux dis- 
tributions positioned on a face of the simplex. To gain 
insight on the pathways involved in the variability of 
our model, we also computed the linear combination of 
fluxes optimized for one of the extreme flux distributions. 
Eight extreme vertices were found for the (Ctrl) and (CN) 
datasets. They were obtained when optimizing the same 
set of linear functions a . In Table 3, the optimal flux distri- 
butions for (Ctrl), (CN) are classified according to the acti- 
vation or inactivation of several pathways in the model. 
It is worth noting that, according to this table, the two 
sets of plausible flux distributions associated with (CN) 
and (Ctrl) have the same topological structure. Indeed, the 
eight extremal behaviors for (CN) and (Ctrl) have clear 
combinatorics: first, NADPH oxidation is set either at 
zero or at its maximal value. Then, a single flux within 
OAA^PYR (R u ), OAA^G3P {R 15 ), G3P^G6P (R 8 ), 
peptide hydrolysis (7? 64 ) is strongly activated whereas the 
three other remaining fluxes are blocked. 

As discussed in a previous paragraph, flux distributions 
with no peptide hydrolysis cannot be considered as rele- 
vant [28,31]. This suggests the extreme flux distributions 
A-F in (CN) and (Ctrl) are not biologically relevant. On 
the contrary, two extreme flux distributions (distributions 
G and H) are consistent with the stoichiometry of the sys- 
tem for the (CN) and (Ctrl) treatments. These two flux 
distributions consist in optimizing peptide synthesis and 
hydrolysis after NADPH oxidation. They corresponded to 
a ratio between peptide hydrolysis and synthesis ^ of 



0.70 in the (Ctrl) treatment and of 0.47 in the (CN) treat- 
ment. These ratios are equal or lower to the maximum 
ratio of 0.67^63 [31]. However, in (CN), treatment peptide 
synthesis was expected to be higher than in Ctrl treat- 
ments since in mammals protein synthesis is reported to 
increase with increasing protein intake (as CN in treat- 
ment) [29]. Curiously, on the contrary, in both flux dis- 
tributions G and H, the total mammary protein synthesis 
decreased for the (CN) treatment when compared to the 
(Ctrl) treatment: V63 equals 430 or 410 mmol/h/half udder 
in (Ctrl) and 283 or 282 mmol/h/half udder in (CN). 

Study of quantitative contributions of precursors (AIO) for 
plausible extreme flux distributions 

As a further investigation to check the relevance of dis- 
tributions G and H in the (CN) and (Ctrl) treatments, 
we studied the quantitative contributions of precursors 
of output nutrients. This study was inspired by the usual 
techniques in the field of nutrition - or any domain con- 
cerned with organ studies. These techniques consist in 
computing yield rates to elucidate how an input nutrient 
may contribute to the composition of an output product, 
for instance to clarify what proportion of glucose, acetate 
or alanine taken up by the mammary gland can be recov- 
ered in the milk components (lactose, fatty acids, protein) 
or oxidized and recovered in C02 released in blood. To 
formalize this issue, we first selected carbon as the compo- 
nent according to which the contributions of precursors 
were to be computed. Then, in order to determine how 
much carbon introduced into the system through a given 
input flux can be recovered in the rate of production of 
an output metabolite, we introduced a precise model for 
the allocation of nutrients in measured outputs (AIO) (see 
Eq. (3)). As shown in Tables 4 and 5, the nutrient allo- 
cation corresponding to both flux distributions G and H 
in the (CN) and (Ctrl) treatments provides evidence that 
glucose is the unique precursor of lactose synthesis. This 
contradicts studies of dairy cows suggesting that glycerol 
and perhaps amino acids contribute to lactose synthesis 
[26], and that fifteen percent of lactose carbon could not 
derive from glucose [27]. 

This precise analysis of the origin of carbon in lactose 
synthesis with AIO suggests that extreme flux distribu- 
tions G and H have to be rejected, although both flux 
distributions are consistent with flux variability criteria. 
We concluded that none of the extreme vertices of the 
set of plausible flux distributions could be considered bio- 
logically relevant with respect to the model and data at 
hand. Notice that the optimization of any linear com- 
bination of metabolic fluxes is either reached by these 
extreme distributions or reached by the infinite number 
of flux distributions lying in a face of the simplex. This 
suggests that the functioning of the mammary gland can- 
not be uniquely modeled by the optimization of any linear 



Table 3 Main properties of the simplex vertices under the assumption of constant ATP-production 
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The qualitative properties of all vertices are shared in (Ctrl) and (CN) treatments. Both correspond to a simplex with height vertices. So are six of each of the (Ctrl) and (CN) simplex vertices. H and G vertices, in the (Ctrl) and 
(CN) treatments, are plausible with respect to the literature. 



Table 4 Origin of carbon mass within outputs for the two optimal flux distributions shown in Table 3 



Origin of carbon mass in outputs for (Ctrl) treatment 
Acetate BHBA Lys Threonine Isoleucine Leucine Valine Histidine Arginine Alanine Glutamine 

Origin of the carbon mass of each output within input (in percentage of total carbon mass of each output) 
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Output 
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Both models have empty flux through the reactions OAA -> PYR (fi, 4 ), OAA-> G3P (R 15 ) and G3P^ G6P (ft 8 ). Model (G) shows strong NADPH oxidation whereas model (H) has zero NADPH oxidation. 
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Table 5 Origin of carbon mass within outputs for the two optimal flux distributions shown in Table 3 



Origin of carbon mass in outputs for (CN) treatment 
Input GLC Glycerol Acetate BHBA Lys Isoleucine Leucine Valine Arginine Glutamine 

Origin of the carbon mass of each output within input 
(in percentage of total carbon mass of each output) 



Output 


Model 


(G) 


(H) 


(G) 




(H) 


(G) 


(H) 


(G) (H) 


(G) 
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(G) (H) 
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0 


0 


0 


0 




0 
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0 




0 


0 


0 


0 
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0 
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0 


0 


0 
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2.9 




1.2 


Peptide 


0 
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0 
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(H) 



Both models have empty flux through the reactions OAA^ PYR (fi, 4 ), OAA-> G3P (R, 5 ) and G3P^ G6P (ft 8 ). Model (G) shows strong NADPH oxidation whereas model (H) has zero NADPH oxidation. 
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combination of metabolic fluxes involved in the current 
model. 

Discriminate treatments despite mammary gland 
flexibility: maxima of AIO on the complete polyhedron of 
flux distributions 

Previous studies suggest that the response of the mam- 
mary gland cannot be modeled uniquely by the optimiza- 
tion of a linear objective function of fluxes. However, 
there are several nonoptimal flux distributions that sat- 
isfy the literature-based information that we have used 
so far b . More generally, many flux distributions compat- 
ible with the additional constraints, including the con- 
dition on carbon precursors for lactose can be shown 
for both (CN) and (Ctrl) treatments. Topological argu- 
ments prove that an infinite set of flux distributions exist. 
Nonetheless, without an idea about the exact shape of 
the space of feasible fluxes (because of the nonlinear 
nature of the condition of carbon precursors), we can- 
not select a plausible point within this space. In other 
words, the available knowledge appears insufficient to 
determine uniquely a flux distribution of nutrients among 
the different branches of the proposed model. 

To understand the functioning of the mammary gland 
and despite this difficulty, we introduced a method to 
estimate the variability of nutrient allocation among path- 
ways on a carbon basis (see Eq. (3)) by computing the 
range (min-max) of AIO coefficients. As these coefficients 
are nonlinear functions of flux variables, computing these 
min and max over the complete space of plausible flux 
distribution requires solving nonlinear optimization prob- 
lems. Our dedicated algorithm detailed in the method 
section scaled properly to the real case that we are study- 
ing 0 . Interestingly, with this approach, we did not favor 
any internal functioning of the system since we parsed all 
objective functions (linear combinations of flux variables), 
in order to have a complete description of the space of 
plausible flux distributions. 

The min-max tables for allocations of nutrients in the 
different pathways provided a clearer view of nutrient 
utilization within the mammary gland. Unlike the func- 
tioning of the extremal distributions (G) and (H) shown in 
Tables 4 and 5, the contribution of carbon from glucose 
to lactose carbon was quite variable within the full set of 
plausible distribution. Indeed, as shown in Tables 6 and 7, 
glucose was the precursor of 85% ( ^ or ^) to 100% 
of the carbon lactose in the (Ctrl) and (CN) treatments. 
Therefore, there exist flux distributions such that glucose 
is not the only precursor of lactose synthesis [26], and 
these distributions are consistent with quantitative esti- 
mations of the ratio of carbon glucose recovered in lactose 
[27]. 

More generally, the intervals of distributions of nutri- 
ents in different pathways were used to compare the 



effects of the (Ctrl) and (CN) treatments. Biologically, in 
comparison to the (Ctrl) treatment, the (CN) treatment 
was characterized by a lower proportion of glucose (on a 
carbon basis) which is oxidized in C02, and a larger ratio 
used for lactose synthesis. This hypothesis can be sus- 
tained since the intervals of distribution of carbon from 
glucose in C02 and lactose almost did not overlap in the 
(Ctrl) and (CN) treatments. More precisely, in the (Ctrl) 
treatment, from ^ = 27.9% to ^ = 39.7% of glu- 
cose carbon was oxidized to produce C02, whereas this 
ratio ranged from 17.5% to 28.9% for the (CN) treatment. 
In addition, in the (Ctrl) treatment, 52.8% to 62.3% of the 
glucose carbon was required to produce lactose carbon, 
whereas 62.2% to 72% of glucose was required during the 
(CN) treatment. 

Altogether, this analysis allowed us to discriminate the 
effects of the different treatments whatever the internal 
functioning of the system may be: the (CN) treatment 
(increase in protein supply to cows) was characterized 
by a lower proportion of glucose oxidized in C02 than 
in (Ctrl). It appears to be a suitable strategy to analyze 
the metabolism flexibility without selecting a precise flux 
distribution or making any assumption on the internal 
metabolic fluxes. 

As a final study, we used an interior point exploration 
method to estimate the AIO variability over the bound- 
ary of the simplex, rather than the complete space that 
was explored previously. As both tables appeared to be 
equal, we concluded that optimized AIO are reached on 
the boundary of the simplex. However, not all extreme 
vertices of the simplex are relevant biologically and they 
do not optimize AIO. This suggests that the flux distri- 
butions that optimize an AIO coefficient are placed at 
the interior of the simplex faces. In the future, it may be 
interesting to study the biological significance with multi- 
objective approaches [34] and check whether these can be 
considered as a "characteristic" point of these faces. 

Impact of long-chain fatty acids oxidation over the model 
predictions 

As a last study, we studied the robustness of our conclu- 
sion with respect to changes in the modeling of input (net 
uptake) of long-chain fatty acids (LCFA), C16 and C18, 
in the tricarboxylic acid (TCA) cycle. In the study of the 
(Ctrl) and (CN) treatments, such an input of LCFA was 
not considered in the system since they are not synthe- 
sized by the mammary gland and we hypothesized that 
they are not oxidized within the mammary gland, based 
on isotope measurements in studies in fed lactating goats 
and in nonruminants [35-37]. This hypothesis was sus- 
tained by the measured carbon balance in the (Ctrl) and 
(CN) datasets [28], which did not require increasing the 
prediction of C02 by introducing any LCFA in the TCA 
cycle (see also Table 8). Nonetheless, in other contexts and 



Table 6 Minimum and maximum utilization of input in each output 



(Ctrl) treatment 

Input Glucose Glycerol Acetate BHBA Lysine Threonine 1 Isoleucine 1 Leucine 1 Valine 1 Histidine 1 Arginine 1 Alanine 1 Glutamine 1 

1422 17,5 1020 336 16,1 1,40 13,1 12,1 12,7 1,38 26,4 9,33 6,10 



Output Minimum and maximum utilisation of Input in each output (in mmol/h/half udder of Carbon) 
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0 




0 


13.4 
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46.6 




0 




0 


37.3 


9.3 




0 




0 




0 




0 




0 




0 






0 






0 




0 


c12 


50.8 




0 
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A local-search algorithm allowed us to compute the minima and maxima of each AIO coefficient for the two treatments (Ctrl), (CN) (in mmol/h/half udder of Carbon). These tables allow discriminating the response of the 
mammary gland to the two treatments without requiring selection of a flux distribution for reactions in the metabolic network. (CN) treatment (protein intake by food) is characterized by a lower proportion of glucose which 
is oxidized in C02 than in (Ctrl). 

(1) Amino acid input corresponded to positive balances between amino acid net uptake and amino acid and utilization in milk protein (i.e. peptide output). 

(2) Amino acid Output corresponded to negative balance between, amino acid net uptake and utilization in milk protein (i.e. peptide output). 



Table 7 Minimum and maximum utilization of input in each output 



(CN) treatment 


Input 


Glucose 


Glycerol 


Acetate 


BHBA 


Lysine 1 


Isoleucine 1 


Leucine 1 


Valine 1 


Arginine 1 


Glutamine 1 




1392 


17.2 


924 


668 


21.5 


21.4 


22.6 


19.3 


26.9 


8.95 
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Minimum and maximum utilisation of input in each output (in mmol/h/half udder of Carbon) 







Min 




Max 


Min 




Max 


Min 


Max 


Min Max 


Min 




Max 


Min 




Max 


Min 




Max 


Min 




Max 


Min 




Max 


Min 




Max 


Glycerol3P 


117 


36.1 




115 


1.8 




1 1.3 


0 


22.8 


0 44.5 


0 




1.6 


0 




2.3 


0.0 




2.0 


0.0 




2.0 


0.0 




2.2 


0.0 




0.9 


Lactose 


1002 


866 




1002 


0 




9.2 


0 


37.9 


0 74.0 


0 




2.6 


0 




3.9 


0.0 




3.3 


0.0 




3.3 


0.0 




3.7 


0.0 




1.5 


C4 


46.4 




0 






0 






0 


46.4 




0 






0 






0 






0 






0 






0 




C6 


33.5 




0 






0 






22.3 


11.2 




0 






0 






0 






0 






0 






0 




C8 


23.0 




0 






0 






17.2 


5.7 




0 






0 






0 






0 






0 






0 




CIO 


64.6 




0 






0 






51.7 


12.9 




0 






0 






0 






0 






0 






0 




C12 


734 




0 






0 






61.2 


12.2 




0 






0 






0 






0 






0 






0 




C14 


250 




0 






0 






215 


35.8 




0 






0 






0 






0 






0 






0 




C16 


343 




0 






0 






300 


42.9 




0 






0 






0 






0 






0 






0 




glycine 2 


6.9 


2.1 




5.5 


0.0 




05 


0.4 


1.3 


0.7 2.6 


0.0 




0.1 


0.0 




0.1 


0.0 




0.1 


0.0 




0.1 


0.0 




0.1 


0.0 




0.1 


Alanine 2 


9.8 


1.4 




9.6 


0.0 




0.9 


0 


2.4 


0 4.7 


0 




0.2 


0 




0.3 


0 




0.2 


0 




0.2 


0 




0.2 


0 




0.1 


glutamate 2 


31.7 


0.8 




8.8 


0.0 




0.4 


6.7 


9.0 


13.0 17.6 


0.4 




0.5 


0.4 




0.7 


0.6 




0.8 


0.3 




0.5 


0.7 




0.9 


0.3 




0.4 


proline 2 


55.0 


1.4 




15.3 


0.0 




0.6 


11.6 


15.6 


22.6 30.5 


0.7 




0.9 


0.8 




1.3 


1.0 




1.4 


0.5 




0.9 


1.2 




1.6 


0.5 




0.7 


aspartate 2 


165 


0.6 




7.0 


0.0 




0.3 


2.7 


4.5 


5.3 8.9 


0.2 




0.3 


0.3 




0.5 


0.2 




0.4 


0.2 




0.4 


0.3 




0.5 


0.1 




0.2 


serine 2 


22.5 


6.9 




22.0 


0.3 




2.2 


0 


4.4 


0 8.5 


0 




0.3 


0 




0.4 


0 




0.4 


0 




0.4 


0 




0.4 


0 




0.2 


C02 


1021 


244 




402 


2.0 




12.5 


180 


227 


351 444 


16.3 




19.8 


14.2 




19.1 


15.8 




20.0 


13.4 




17.4 


14.8 




19.2 


5.9 




7.8 



A local-search algorithm allowed us to compute the minima and maxima of each AIO coefficient for the two treatments (Ctrl), (CN) (in mmol/h/half udder of Carbon). These tables allow discriminating the response of the 
mammary gland to the two treatments without requiring selection of a flux distribution for reactions in the metabolic network. (CN) treatment (protein intake by food) is characterized by a lower proportion of glucose which 
is oxidized in C02 than in (Ctrl). 

(1) Amino acid input corresponded to positive balances between amino acid net uptake and amino acid and utilization in milk protein (i.e. peptide output). 

(2) Amino acid Output corresponded to negative balance between, amino acid net uptake and utilization in milk protein (i.e. peptide output). 
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Table 8 Effect of long-chain fatty acids (LCFA) oxidization in the triacarboxylic acid (TCA) cycle over model analysis 



Main characteristics of models with different ratios of LCFA oxidized in TCA 



Model 






C02 prediction 


ATP 


ATP = 1 250 mmol/h/half udder 














optimum 


CXllcmc TIUX CHbllllJUllun Ibee 1 dule 




Variability of AIO coefficients 


Ratio of 
long-chain FA 
oxidated in TCA 


Dataset 


Predicted Ratio of 

C02 predicted C02 

and measured C02 




Total number Nonplausible 
flux values 


Nonplausible 
AIO 


UlULUac Lai UUI 1 

required to 
produce lactose 


uiulujc Laiuuii 

oxidized 




(HB) 


1546 


Non available 


6628 


Nonrelevant hypothesis [20-22] 








0% 


(CN) 


1021 


99% 


2045 


8 6 


2 


[62.2 ; 72.0] 


[17.5; 28.9] 




(Ctrl) 


1126 


121% 


3081 


8 6 


2 


[52.8 ; 62.3] 


[27.8 ; 39.7] 




(HB) 


1640 


Non available 


7395 


13 13 


0 


[47.7 ; 62.4] 


[28.8; 45.9] 


10% 


(CN) 


1107 


107% 


2739 


8 6 


2 


[59.4 ; 72.0] 


[17.5; 32.0] 




(Ctrl) 


1202 


129% 


3701 


8 6 


2 


[5 1.6; 62.3] 


[27.9 ; 41 .2] 




(HB) 


1756 


Non available 


8336 


13 13 


0 


[46.9 ; 62.4] 


[29.0; 47.0] 


20% 


(CN) 


1214 


1 1 8% 


3607 


8 6 


2 


[57.3 ; 72.0] 


[17.5; 34.6] 




(Ctrl) 


1298 


140% 


4476 


Nonrelevant hypothesis: predicted C02 is not compatible with measured C02 [28] 








(HB) 


1826 


Non available 


8396 


13 13 


0 


[46.5; 62.4] 


[29.0; 47.5] 


25% 


(CN) 


1279 


1 24% 


4128 


8 6 


2 


[56.3 ; 72.0] 


[17.6; 35.8] 




(Ctrl) 


1355 


146% 


4938 


Nonrelevant hypothesis: predicted C02 is not compatible with measured C02 [28] 







ATP balance 


Extreme distributions are 


For plausible ratios of long-chain 


is too high 


not biologically relevant 


FA in TCA, (CN) treatment is 






characterized by a lower proportion 






of glucose (on a carbon basis) 






which is oxidized in C02, and 






a larger ratio used for lactose synthesis. 



To study the impact of the variability of FA oxidation, a ratio of long-chain FA (10%, 20%, 25%) was introduced in the TCA cycle. For datasets (CN), (Ctrl) and (HB) which were compatible with a given ratio of LCFA oxidation, 
extreme flux distributions and AIO coefficients variability were studied. Our main conclusions are robust to the introduction of LCFA oxidation (Table 8). Interestingly, as shown in Table 9, the structure of the simplex 
generated by the (HB) diet is more complicated than the (Ctrl) and (CN) treatments, with 13 vertices for all hypotheses. This is due to the fact that R\s is highly constrained by measurements, so that this flux cannot vanish 
when R 14 is optimized or when R K is maximized. All the vertices for the (HB)-simplex contradict knowledge-based literature. 
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models, LCFA oxidation was either measured in starva- 
tion condition [38] or introduced in order balance carbon 
consumption and production with respect to C02 pre- 
diction [18,20-22]. This LCFA oxidation hypothesis was 
for instance retained in models studying the (HB) dataset 
[20,21]. 

From this literature study, it appears that LCFA oxida- 
tion may depend on both the environmental and experi- 
mental contexts, and no single model can be favored yet. 
To study the impact of LCFA oxidation, we successively 
introduced a ratio of LCFA in the TCA cycle (10%, 20%, 
25%), and assumed that 50% of C16 output are synthesized 
within the mammary gland [21] (see Figure 3). Then we 
performed a complete study of the (Ctrl), (CN) and (HB) 
datasets for these hypotheses (see Tables 8 and 9). 

For the (Ctrl) treatment, by comparing the predicted 
and measured C02 quantities, we concluded that the 
hypotheses of 20% and 25% of FA oxidized in the TCA 
cycle were not in agreement with our experimental data: 
the increase of C02 prediction was too large both when 
compared to measured C02 (see Table 8) and when con- 
fronted to the measured increase of 9% of C02 deriving 
from long LCFA in lactating goats in an extreme star- 
vation condition [38]. Similarly, the hypothesis of 0% of 
LCFA oxidation was not consistent with the (HB) treat- 
ment [22]. 

For all remaining compatible pairs of model and 
datasets, we first studied the ATP maximization hypothe- 
sis. In all cases, our results, shown in Table 8, suggest that 
ATP maximization is not biologically relevant: ATP bal- 
ance was even larger than the ATP balance of the models 
where no LCFA was introduced in the TCA cycle (0%). 



Then we enumerated extreme flux distributions and 
studied their biological relevance. As shown in Table 9, 
for each ratio of LCFA in the TCA cycle, the struc- 
ture of the simplex of plausible flux distributions for the 
(Ctrl) and (CN) datasets was similar to that shown in 
Table 3 (0%). Using similar arguments, no extreme dis- 
tribution could be considered as biologically relevant. 
When studying the space of solution associated with the 
(HB) dataset, a more complex structure appeared, based 
on 13 extreme flux distributions instead of 8. Explicit 
(although not unique) linear functions that can be opti- 
mized to obtain these extreme distributions are detailed in 
Table 9 d . More precisely, for the (HB) treatment, we still 
obtain a division according to the activation of the four 
fluxes OAA^PYR (i? 14 ), OAA^G3P (R 15 ), G3P^G6P 
{Rs), peptide hydrolysis (i?64)- Nonetheless, an intricate 
phenomenon appears. Indeed, measurements imply that 
OAA^G3P (Ris) is very constrained. It cannot vanish at 
the same time as OAA->PYR (7? 14 ) or when NADPH oxi- 
dation (i?i9) is at its maximal value. Therefore, for several 
optimal conditions related to NADPH oxidation (i?i9), 
OAA^PYR(i?i 4 ), G3P^G6P (R 8 ), peptide hydroly- 
sis (i?64)> we need to decide whether OAA^-G3P (R15) is 
slightly non-zero, or whether OAA— >G3P is assumed to 
vanish although the optimal is not exactly reached for the 
initial flux considered. Despite this difference of struc- 
ture between the (HB) dataset and the (Ctrl) and (CN) 
datasets, our analysis suggested that extreme distributions 
were not biologically relevant, independently from the 
ratio of LCFA oxidation or the dataset under study. 

Finally, comparing the variability of the AIO coefficients 
for the (Ctrl) and (CN) treatments when oxidizing 10% of 




T NADHc / 

NH3 I, , R_24 jT / 



Figure 3 Modeling long-chain fatty acids oxidation (LCFA) in tricarboxylic acid (TCA) cycle. Introducing LCFA (C(16:0), C(18:0)) oxidation in 
the TCA cycle may be required to consistently model the response to several treatments, such as the (HB) dataset (Table 1). In this case, the model 
shown in Figure 2 is extended by introducing inputs of C(16:0) (ftm) and C(18:0) {Ru2), and output of C(18:0) (Rus)- C(16:0) oxidation (R20) and 
C(1 8:0) oxidation (R21 ) are modified accordingly. 



Table 9 Effect of long-chain fatty acids (LCFA) oxidization in the triacarboxylic acid (TCA) cycle over model analysis 



Main properties of the simplex vertices, assuming constant ATP-production, 
with different ratios of LCFA oxidized in TCA 





% of FA 
oxidated 
in TCA 


name 


Fy^mnlp nf 

LAQI 1 1 1-/ 1 d \Jl 

maximized 
function 






fnmhinatnrir«; 

of pathways 






Validation 












R 19 
























NADPH 
oxidation 


OAA^G3P 


OAA-* PYR 


G3P^G6P 


Peptide 
hydrolysis 


Peptide 
synthesis 


Pyr ->- OAA 




(Ctrl) 


0% 
10% 








1831 
2451 








125 


1835 
2455 




(CN) 


0% 
10% 
20% 
25% 


B 


V15-V19 


0 


795 
1489 
2357 
2878 


0 


0 


0 


150 


803 
1497 

2365 
2886 


Non 

relevant 

flux 

values 

for 




10% 








6173 










6145 


(HB) 


20% 
25% 








7115 
7675 








150 


7086 
7646 




(Ctrl) 


0% 
10% 










1831 
2451 






125 


1835 
2455 




(CN) 


0% 
10% 
20% 
25% 


F 


1/14-1/, g 


0 


0 


795 
1489 
2357 
2878 


0 


0 


150 


803 
1497 

2365 
2886 


Non 

relevant 

flux 

values 

for 

Rl3,ft64 




10% 










6173 








6145 


(HB) 


20% 
25% 










7115 
7675 






150 


7086 
7646 




(Ctrl) 


0% 
10% 












3662 
4902 




125 


4 




(CN) 


0% 
10% 
20% 
25% 


D 


V 8 -V 19 


0 


0 


0 


1590 
2978 
4714 
5756 


0 


150 


8 


Non 

relevant 

flux 




10% 

20% 


D1 


V 8 -V,9-V 14 




29 


0 


12289 
14172 








values 
for 

ft8,/?64 


(HB) 


25% 






0 






15292 


0 


150 


0 





Table 9 Effect of long-chain fatty acids (LCFA) oxidization in the triacarboxylic acid (TCA) cycle over model analysis (Continued) 





10% 












12289 












20% 


D2 


V8-V1 9-^15 




0 


29 


14172 












25% 












15292 










(Ctrl) 


0% 
10% 














305 
409 


430 
533 


4 


Glucose is 
the unique 
precursor of 
lactose 
synthesis 
(AIO) 


(CN) 


0% 
10% 
20% 
25% 


H 


1/64 -V: 9 


0 


0 


0 


0 


133 
248 
393 
480 


283 
398 
543 
630 


8 


(HB) 


10% 

20% 
25% 


H1 


V(, 4 -Vig-Vi4 


0 


29 


0 


0 


1024 
1181 
1274 


1174 
1331 
1424 


0 


Non 
relevant 
flux 
values 
for R 63 


10% 














1024 


1174 






20% 


H2 


V64-V] 9-^15 




0 


29 




1181 


1331 








25% 














1274 


1424 






(Ctrl) 


0% 
10% 
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694 


1714 

2330 








125 


1718 
2334 




(CN) 


0% 

10% 

20% 


A 


V15+V19 


22 


791 

1485 

2353 


0 


0 


0 


150 


799 

1493 

2361 


Non 

relevant 

flux 

values 

for 


Extreme flux 
distributions 

within the set of (HB) 

plausible 

solutions 


10% 

20% 
25% 






1216 


5961 
6902 
7462 








150 


5932 
6873 
7433 
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0% 
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694 




1714 

2330 
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0% 
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E 
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22 
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Table 9 Effect of long-chain fatty acids (LCFA) oxidization in the triacarboxylic acid (TCA) cycle over model analysis (Continued) 



(Ctrl) 



(CN) 



(HB) 



(Ctrl) 



(CN) 



(HB) 



10% 

20% 

25% 

0% 

10% 

0% 

10% 

15% 

20% 

10% 

20% 

25% 

10% 

20% 

25% 

0% 

10% 

0% 

10% 

20% 

25% 

10% 

20% 

25% 

10% 

20% 

25% 



E2 1/14+V19-50V15 



CI 



Gl 



Vg+Vig 



C2 1/8+V19-50V15 



V64+V\9 



V64+V\9 



G2 1/64+Vl 9-501/1 5 



Litterature-based upperbounds for fluxes 



946 



694 



22 



1216 



946 



694 



22 



1216 



946 



32 



32 



6008 
6949 
7509 



29 



3428 

4659 

1583 

2971 

4706 

5749 

11858 

13840 

14961 

11958 

13840 

14961 



29 



mmol/ 

h/half udder [33] 



< 591 
[28,31] 



125 



150 



286 

388 

132 

248 

392 

479 

989 

1145 

1238 

996 

1146 

1247 

Non-zero 
whole body 
protein synthesis [29] 



450- 



410 

513 

282 

398 

542 

629 

1139 

1295 

1388 

1146 

1296 

1397 



Lower than 
udder [33] 



5979 
6920 
7481 



< 266 mmol/h/half 



Non 

relevant 

flux 

values 

for 



Glucose is 
the unique 
precursor of 
lactose 
synthesis 
(AIO) 



Non 

relevant 

flux 

values 

for 

f?63,R64 



To study the impact of the variability of FA oxidation, a ratio of long-chain FA (10%, 20%, 25%) was introduced in the TCA cycle. For datasets (CN), (Ctrl) and (HB) which were compatible with a given ratio of LCFA oxidation, 
extreme flux distributions and AIO coefficients variability were studied. Our main conclusions are robust to the introduction of LCFA oxidation (Table 8). Interestingly, as shown in (Table 9), the structure of the simplex 
generated by the (HB) diet is more complicated than the (Ctrl) and (CN) treatments, with 13 vertices for all hypotheses. This is due to the fact that ft 15 is highly constrained by measurements, so that this flux cannot vanish 
when Ri4 is optimized or when S19 is maximized. All the vertices for the (HB)-simplex contradict knowledge-based literature. 



J > 
3 S" 

TJ Q- 

O 

r> 

b~ ns 

§ ° 



O NJ 

Ln O 

0 T* 

ID i» 

£» CO 

00 00 



Abdou-Arbi etal. BMC Systems Biology 2014, 8:8 
http://www.biomedcentral.eom/1752-0509/8/8 



Page 19 of 26 



LCFA in the TCA cycle still suggested that the (CN) treat- 
ment is characterized by a lower proportion of glucose (on 
a carbon basis) which is oxidized in C02, and a larger ratio 
used for lactose synthesis. 

Altogether, this study suggests that the main character- 
istics of the (Ctrl) and the (CN) treatments are robust and 
could be elucidated despite lacking information on the 
precise internal behavior of LCFA oxidation. 

Discussions 

Analyzing the distribution of nutrients in a metabolic 
network to study the flexibility of a metabolism at the 
organ level 

An important challenge in applicational fields of 
metabolism studies at the organ level is to understand 
how the components of inputs are transformed into some 
expected outputs, under some assumptions about the 
functioning of a system. To that end, great use is made of 
comparisons between yield rates describing the allocation 
of input nutrients within the set of outputs. Nonetheless, 
to allow a precise comparison of nutrients, these studies 
require insights on the distribution of matter components 
across the output of each reaction involved in the system. 
Such information may be provided by several experi- 
mental marking techniques (C13 MFA) that make use of 
carbon isotopes [39]. Nonetheless, their application to 
mammals is challenging [40]. 

To elucidate how input nutrients are allocated among 
the output nutrients of a metabolic system despite exper- 
imental limitations, we have introduced novel methods 
which refine the flux balance analysis of a metabolic sys- 
tem related to an organ of a large animal. Our method 
can be seen as an extension of Flux Variability Analysis 
[11], where the emphasize is put on efficiency or yield rate 
variability rather than on flux variability. 

As an example of application, we have studied the mam- 
mary metabolism in ruminants (dairy cows) [15,23,24]. 
Compared to the conventional models available, our sto- 
ichiometric model describes very precisely the variations 
of energy consumptions as ATP, allowing us to inves- 
tigate some optimization hypotheses related to energy 
variations. 

As a methodological innovation, we introduced a 
method to estimate the Allocation of Inputs in Outputs 
(AIO), that is, the ratio of transformation of each input 
nutrient into outputs, provided that we are given a flux 
distribution balancing the production and consumption 
of intermediary metabolites. Two AIO computations are 
possible: one may fix a linear combination of flux variable 
to optimize, leading to an extreme flux distribution (i.e. 
extreme pathways) whose allocations (AIO) are precisely 
computed with our method. Alternatively, we can rea- 
son with regard to the possibly infinite complete convex 
space of plausible flux distributions, in the spirit of FVA 



analysis, without solving the set of equations. In this case, 
our method estimates the minima and maxima of each 
allocation (AIO) coefficient within the complete space. 

The latter point is where our main algorithmic inno- 
vation lies. Indeed, AIO coefficients are not linear with 
respect to the flux variables v. Therefore, they have no rea- 
son to be extremal for flux distributions corresponding to 
vertices, edges or faces of the simplex, unless the gradi- 
ent of every component of the AIO is a nonzero function. 
Therefore, we have introduced two algorithms allowing us 
to compute the extremal values of AIO coefficients when 
the flux space is bounded. The fastest algorithm requires 
the exhibition of an analytic expression of the AIO matrix. 
If this is not possible, a local-search algorithm can be 
used. These algorithms can be used to check whether the 
extremal values are reached for flux distributions lying on 
the boundary of the simplex of plausible flux distributions. 

This approach applied to a stoichiometry model taking 
ATP variations into account permitted a better under- 
standing of ruminant mammary gland metabolism in 
comparison to previous studies based on similar mod- 
els and datasets [22,24]. For instance, we were able to 
characterize the differences in the effects of two treat- 
ments, such as the quantitative range of the proportion 
of glucose which is oxidized in C02, or used for lactose 
synthesis. This information was derived from the carbon 
composition of each metabolite, expert knowledge which 
is easily accessible and is therefore compatible with the 
experimental limitations regarding mammals. 

Optimization strategies within tissue or organs 

Our study allowed us to revisit optimization-based 
hypotheses on the functioning of the mammary 
metabolism. We have provided evidence that flux distri- 
butions corresponding to an optimal production of energy 
(ATP) cannot describe appropriately the metabolism of 
the mammary gland, as would be the case for bacterial 
metabolisms, which tend to optimize biomass-related 
functions [5,41]. The main reason is that the ATP balance 
is too large and variable among the responses to different 
treatments, confirming previous observations [31]. As a 
consequence, the underlying hypotheses used to drive 
previous studies of mammary metabolism have to be 
carefully reformulated [15,21,23]. 

As an alternative, we have hypothesized that ATP bal- 
ance remains nonoptimal and almost constant in response 
to several treatments, and we have introduced in the 
model a recent estimation for this quantity [31]. Our goal 
was to check whether the observed responses to the sys- 
tem could be explained by the optimization of a linear 
combination of fluxes, that is, an extremal vertex of the 
simplex of plausible flux distributions. This hypothesis 
was rejected, first because it led to non-realistic orders of 
magnitude for some fluxes such as peptide hydrolysis, and 
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second because the precursors of some components such 
as glucose were not biologically relevant [26,27]. Notably, 
it was necessary to trace back the quantitative origin of the 
nutrients to complete this analysis and reject all extreme 
flux distributions, illustrating the advantage of the AIO 
approach. 

Therefore, no optimization of a linear combination of 
flux variables could be found to uniquely describe the 
metabolism of the mammary gland, a multicellular organ, 
as an extreme flux distribution of our model. To overcome 
this limitation of FBA-inspired analyses at the mammalian 
tissue or organ level, we studied the variability of AIO 
coefficients by introducing the computation of min-max 
ranges of AIO. This led to a general overview of the effects 
of treatments without precluding any steady state internal 
behavior. 

Benefits of studying the variability of the allocation of 
input in output (AIO) and future model refinements 

The study of the variability of AIO, that is, intervals of 
allocation of nutrients, made it possible to distinguish the 
metabolism of the mammary gland in two different nutri- 
tional conditions corresponding to an increase in protein 
supply (Ctrl and CN). Notably, the intervals of allocation 
of glucose in C02 and lactose were different in the two 
nutritional conditions, reflecting the mammary gland's 
metabolic flexibility. Nonetheless, several improvements 
can be made to the model. 

First, the simplex of plausible distributions could be 
reduced by introducing knowledge currently available 
on the kinetic bounds for enzymatic activities, intro- 
duced in previous numerical models [21]. This may for 
instance prevent some non-constrained behaviors shown 
in Tables 6 and 7, such as the possible nonzero contri- 
bution of acetate and /J-hydroxybutyrate (BHBA) carbon 
to lactose. The main issue in this area will be to aggre- 
gate the enzymatic knowledge to fit with the format of 
the reactions in the model, such as the reaction which 
transforms G3P to PYR under the regulation of five differ- 
ent enzymes, together with NAD+ and ATP availability. 
Therefore, proposing relevant maximal values requires 
the use of advanced methods for model reductions. 

A second improvement of the model is dependent on 
the production of additional observations to clarify the 
set of possible behaviors of the system. Particularly, the 
contribution of acetate and /i-hydroxybutyrate (BHBA) 
to lactose shown in Tables 6 and 7 remains question- 
able in ruminants since it suggests that the carbon of 
acetyl-COA could leave the TCA cycle of the mitochon- 
dria through the OAA G3P pathway. This could be 
explained by the fact that the carbon allocation in the 
model did not take account of the carbon positions. How- 
ever, one isotope model and one study of gluconeogenesis 
in man suggested that labeled carbon within glucose came 



from [2 - 14C]-acetate (through 14C-acetylCoA in TCA) 
[42,43]. It remains challenging to obtain such precise 
information for large ruminants such as dairy cows. Nev- 
ertheless, the model could be improved by modeling with 
additional constraints the few items of knowledge we have 
on the flux distribution of positioned carbons. 

The last improvement of the model consists of includ- 
ing constraints that do not correspond to rule-based 
metabolic equations. More precisely, in several numeri- 
cal models, the effects of external fluxes were introduced 
to predict the response of the mammary gland in a con- 
sistent way [14,17,18]. Among external fluxes, we noticed 
the regulatory effect of long-chain fatty acids on the activ- 
ity of enzymes involved in fatty acid synthesis, although 
the mechanisms are not clearly understood [44]. Although 
some extensions of the FBA formalism to dynamic reg- 
ulation exist [45], the data availability was insufficient to 
apply the extensions in the model. Therefore, existing reg- 
ulations of enzyme activities cannot be included directly 
in the stoichiometric model. To overcome this limitation, 
we plan to investigate the effects of external fluxes by 
introducing additional constraints in the model, to check 
whether the nutrient output in the milk can be predicted 
from the nutrient input. 

Conclusion 

We have introduced a method in the framework of flux- 
balance analysis of a metabolic network. As a main nov- 
elty, our approach allows studying the variability of effi- 
ciencies (or yield rates) of a metabolic model provided 
with input-output measurements. More precisely, our 
approach allows a quantitative estimation of the minimum 
and maximum proportions of the carbon quantity of each 
input nutrient which is recovered in each output compo- 
nent of the system. The main innovation is to propose a 
method which does not require determining the quantita- 
tive distribution of nutrients between the branches of the 
system. To that end, we have performed a parsing of the 
space of flux distributions which are compatible with both 
the model stoichiometry and input-output measurements. 

This method was applied to study the response of the 
mammary gland to several treatments. It allowed us to 
distinguish two different metabolic responses of the sys- 
tem, corresponding to two nutritional situations and accu- 
rately reflecting metabolic flexibility. Overall, our method 
appears to be configured to study the variability of the 
yield rates of a metabolic system at the multicellular or 
organ level without making any hypothesis on the internal 
behavior of the system. 

Method 

Flux modeling based on Flux Balance Analysis 

Metabolic models are described according to the generic 
framework of Flux Balance Analysis [6,25,46]. Metabolites 
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split into three non-intersecting sets. First, input metabo- 
lites are gathered in set I. The rate vector of all input 
fluxes (with in the form /") is denoted by vj € W. Sec- 
ond, output metabolites are denoted by O. We denote by 
vo € R ? the rate vector of the corresponding output fluxes 
(in the form "O -*■"). Finally, intermediary metabolites (or 
pivots) are denoted by V, and its cardinality is denoted by 
n. The complete set of metabolites isA4=2UVUO = 
{mr . ..m p+n+q ). 

Let 1Z = {ri . . . r t } be the set of t reactions which pro- 
duce some metabolites while degrading others. The rate 
vector of these reactions is denoted by v e R l . A reac- 
tion r) has the form £ mGluT > s m,j m — ► Hm'&vuo s 'm,j m ' 
where s m j and s' m j are input and output stoichiometric 
coefficients. The metabolic system is described by the so- 
called stoichiometric matrix M = (aij)i< p+tl+q j< r , where 

a hi = s ' mi ,j ~ s m,j- 
As a usual assumption, intermediary metabolites cannot 

accumulate in the cell: the consumption and production 

flux rates of intermediary metabolites are balanced. A 

linear constraint is derived on the rate vector v e R f . 

Additional biological knowledge about the distribution 

of fluxes within the system is modeled by a matrix B. 

Overall, plausible flux distributions v satisfy the following 

equations: 

/-vA 

Ml v I = 0 [Stoichiometry constraints]; 
\ v 0 / 

Bv = 0 [Biological knowledge constraints]. 

(1) 

Assuming that reactions are irreversible and fluxes have 
physical upper-bounds, the set of solutions to Eq.(l) is 
a convex polyhedron, called a simplex, of plausible flux 
distributions. The simplex can be described by means of 
its vertices, edges or faces [6,25,46], which are related to 
the optimization of objective functions and have shown 
their biological significance in several contexts [47]. In 
the example described below, vertices of the simplex (i.e. 
extreme pathways) were computed with standard linear- 
based methods in the bound case, and with probabilistic 
methods in the unbound case. As a refinement of these 
methods, we used a random sampling of the set of possible 
linear objective functions to obtain a complete descrip- 
tion of the set of extreme pathways. This approach was 
efficient since the dimension of the space of feasible dis- 
tribution is quite small, because it is strongly constrained 
by the input-output vectors vj and vq. 

Modeling the quantitative contribution of input 
metabolites to output nutrients: AIO 

In studies at the organ level such as in nutrition, the 
choice of a plausible flux distribution v aims to elucidate 



how an input nutrient may contribute to the composi- 
tion of an output product. To formalize this issue, we 
introduce a component according to which the allocation 
of input nutrients will be computed. For dietary applica- 
tions, carbon is a relevant element since it appears in the 
composition of all metabolites in the system. Nitrogen is 
less generic but more relevant when specifically studying 
nitrogen metabolism. Our issue therefore reads as follows: 
how much carbon introduced into the system through a 
given input flux can be recovered in the rate of production 
of an output metabolite? 

We introduce the following formalism. Let c € W +n+q 
be a vector describing the component composition of 
metabolites (its carbon composition for instance). Assume 
that the stoichiometry of each reaction in the system 
satisfies a matter-invariant property with respect to the 
component, that is: J^rmeM a i,j c ( m i) = 0 f° r every reac- 
tion rj € 1Z. Of the total component mass provided as 
a substrate to a reaction r, only a part contributes to the 
composition of a given product m. Let In r {m) denote this 
ratio. Assuming that the system follows a proportional 
matter distribution, a consequence of the mass-invariant 
property is that XimeA-l = 1 f° r a U r € 1Z. 

Let v be a plausible flux distribution and a, solution to 
Eq.(l). Let m be an input metabolite, for example glu- 
cose. The rate of the flux of component mass brought 
into the system by m equals C(m)vj(m). We denote by 
Xo[v, m] 6 the vector of proportions of input compo- 
nent fluxes recovered in each output flux (for instance, the 
proportion of carbon from the input glucose appearing in 
the composition of the C02 in blood). The coefficients of 
xo[v,m\ are called the allocation of the metabolite m in 
output nutrients. 

To determine these ratios, we introduce xp[v,m] e W 
the vector of ratios of input component fluxes which 
are required to produce each intermediary metabolite 
m before its breakdown into other metabolites. The 
law of matter conservation and the assumption that 
intermediary metabolites do not accumulate put sev- 
eral constraints on xp[v,m] and xo[v,m], which are 
detailed in Tables 8 and 9. We deduce that the vec- 
tor of allocations xo[v,m\ is a solution to the equation 

D 2 [v]x ( / n) = -^iM^^'m]) for the input param- 
eter vector xf 1 ^ = (0, . . . , 0, C(m)vj(m), 0, . . . , 0), where 
D\ and D 2 are defined in Table 10. To elucidate whether 
this system of matrices determines uniquely xol v, m], we 
noticed that D\ is a square matrix such that all diagonal 
coefficients are equal to 1 and the others are nonpositive. 
This family of matrices is named the M-matrix in [48], 
and it has been shown that such a matrix is invertible. We 
deduce that D\ is invertible so that the vectors of allo- 
cations Xo[v,m] and Xp[v,m] are determined uniquely. 
With these formulas, when the flux distribution v has been 
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Table 10 Modeling the quantitative allocation of input nutrients in output products 

ln'(m k ) 



v akf{m i w< s ifow > 0 



0 



otherwise. 



Ratio of the flux of component c provided as substrate to a 
reaction rj recovered in the composition of the product m k . It 
is the sum of individual substrate contributions. 



F ( m i) - J2r n e'R. l a Ln >0 a i,ri v 1 ~ Hr, e"R,a, ( , <0 Kl K 


Total metabolite rate involved in the production of an interme- 
diary metabolite m,-, before its degradation by other reactions 


dkj 


1 \fk=i 
- £ n^*ln r i(m k ) otherwise. 


Ratio of a product flux (m^ 6 PUO) on the production of 
m, e I U P U 0 


Ddv] 
D 2 [v] 


(dkj)p<kj<p+n+q 
(dkj)p<k<p+n+q,l </<p 


Linear transformation of matter components contained in 
intermediary or output metabolites 

Linear transformation of matter component contained in input 
metabolites 


x[ v, m] = 


( *' m) l 
x P [x,m] 

v x 0 [x,m] / 


Rates of fluxes of component brought by the m-input flux 
appearing in the composition of each metabolites. x' m ' = 
(0 0,C(m)v/(m),0 0). 


D 2 [v]xf } 


/ xp[x,m] \ 

-D,[v] 

yx 0 [x,m] y 


Constraints on component fluxes deduced from the matter- 
invariance law, derived from Eq.(2) below. 



We are given a stoichiometric matrix A and a vector c e R p+n+<1 which describes the component composition of all metabolites. If v is a fixed flux distribution which is 
compatible with the stoichiometry of the system, AIO[ v] is a matrix whose (/",/) input describes the proportion of component quantity contained in m, which is 
recovered in the flux of the output m. 



x[v,m](m k )= ^ 



ln'i(m k ) 



F(m f ) 



x[ v, m] (mi) 



(2) 




proportion of component ■ ' rate of the flux of 

substrate of/j Proportion of 



flux brought to the 
composition ofm^ 
through/) 



rrij — component 
consumed by/) 



component 
appearing in 
the composition ofm, 



total flux of component provided as substrate tor, 



chosen, any algorithm or method to inverse matrix D\ can 
be used. 

Overall, the complete table showing the distribution of 
nutrient inputs in nutrient outputs is defined to be the q x 
p matrix AIO [ v], the columns of which are Xo[ v, m{\, . . . , 
xq\ v, m p ], as follows: 

AIO[ v] = xq[v, m] . 



\xo[v,m\ ) 



' C(m\)vi(m\) 



C(m p )vi(m p ) / 



(3) 



With these formulas at hand, as soon as the flux distribu- 
tion v has be chosen, any algorithm or method to inverse 
the matrix D\ can be used. Altogether, the computation of 



a complete AIO table takes 0((n+p+q) 3 + np) operations 
for a given flux distribution v. 

Computing extrema of the AIO coefficients among the 
complete polyhedron of plausible flux distributions: 
solving nonlinear optimization problems 

Computing the extrema of the AIO coefficients among the 
convex polyhedron of plausible flux distributions v, leads 
to an optimization problem in a nonlinear context. The 
problem rewrites in an optimization scheme, leading to 
2q x p optimization problems of the form 



Vj < q,Vk < p, (Max/Min)imize 

/-V 

AIO[v] [i,k]; v e K", M I v 



= 0, Bv = 0, v > 0 



(4) 
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These problems are nonlinear programming (NLP) prob- 
lems [49] since they aim to optimize a nonlinear objec- 
tive function over a possibly nonbounded simplex search 
space. 

When both the simplex search space was bounded and 
an analytical expression for AIO [ V] as a function of all the 
free variables of the eq.(l) solution space was producible 
with formal algebra software [50], the Matlab fmincon 
function was used to solve the NLP problem. Among 
all available optimization routines, the best performance 
was achieved by the interior point algorithm, which has 
a polynomial complexity for these problems [51]. Other 
optimization routines were used to confirm our results, 
but the computation time was much longer. 

When an analytical expression was not available, we 
implemented an alternative strategy making use of some 
local search routines together with the fact that AIO(v) 
can be computed for any v in a cubic time. First, we 
compute the Chebyshev ball B(xo,r) of the simplex [52]. 
Any point x of the w-dimensional ball is uniquely iden- 
tified by a vector (0i,...,0„) with 0; €[0,2tt[ for 1 < 
i < n and 9 n e [ 0, it [. The intersection P(x) of the half 
line [#o>*) with the simplex of plausible flux distribu- 
tions is unique if it exists. Conversely, for any point y 
on the boundary of the simplex, there exists a unique 
x € B{xq, r) such that y = P(x). By extension, we say 
that (0i, . . . , 0„) is a parameterization of the boundary of 
the simplex. This parameterization is used to determine a 
discretization of the simplex boundary. Given an integer 
p > 0, every 0,- is assumed to be of the form 2jzj/p, for 
j € {0, . . . ,p — 1}, whereas 9 n rewrites as Ttj/p. The finite 
neighborhood of any point y on the boundary, denoted 
by Af(y) is obtained by modifying each single coordi- 
nate of the parameterization {6\, . . ., 2nj/p, ...,9 n ) of y: 
for instance, among the maximal 2n neighbors, we find 
the point with paramerization (0i, . . . ,2n(j— l)/p, . . . , 9„) 
and (0i, . . . , 2n(j + l)/p,. . ., 9„). We can finally apply a 
local search algorithm based on a dichotomy principle to 
improve the discretization level and reach the solution to 
the optimization problem on the simplex boundary. 

The same algorithm is extended to the full space by 
introducing a supplementary coordinate standing for the 
distance between P(x) and the center Xo of the Chebyshev 
ball. This makes it possible to check whether the optimum 
identified previously lies on the boundary of the simplex. 

Analysis workflow 

A web application dedicated to the computation of AIO 
for metabolic networks is freely available [53]. It enables 
us to build a stoichiometric model, define its inputs and 
outputs, solve the associated FBA problem and finally 
compute the AIO of the selected flux distribution. Offline 
tools enabling the computation of the extrema of the AIO 
are also provided in a companion web page 6 . 



The complete workflow of analysis requires the online 
webpage together with the use of several environments. 
First, php is required for the computation of extremal ver- 
tices of the convex polyhedron of plausible flux distribu- 
tions, including the case when this space is not bounded. 
Second, SAGE is used to performed formal algebra tasks 
involved in the computation of AIO matrices. Finally, 
matlab is needed to compute the extrema of the AIO 
coefficients. The main functionalities of the workflow are 
depicted in Figure 2. 

Mammary gland stoichiometry model 

As the basis of the ruminant mammary gland metabolism 
that is studied in this paper, we used a generic model of 
mammary metabolism [15,23]. This generic model con- 
tained 54 reactions involving six intermediary metabo- 
lites, called carbon-chain pivots, at cross-over points 
between metabolic pathways, chosen to describe the main 
possible conversions between metabolites: oxaloacetate 
(OAA), a-ketoglutarate (a-KG), pyruvate (PYR), acetyl 
coenzyme A (AcoA), glucose (GLC) and serine (SER pivot). 
The reactions also take into account the variations in ATP, 
four cofactors (NADHc, NADHm, NADPH and FADH2) 
and other metabolites (02, C02 and NHS). 

The model was first extended and detailed in order 
to include the reactions included in other models of 
mammary metabolism [15,21,22,24]. Reactions for lac- 
tose synthesis, milk protein, fatty acids and glycerol-3P of 
triglycerides were included in the model, as well as four 
new pivots. Mitochondrial (mAcoA) and cytosolic (cAcoA) 
acetyl-CoA were used to distinguish between utilization 
of cACoA for fatty acid (FA) synthesis and utilization of 
mACoA in the tricarboxylic acid (TCA) cycle. In addition, 
glucose-6-phosphate (G6P) and triose-phosphate (G3P) 
were introduced to account more specifically for the 
transformation of glucose through the pentose phosphate 
pathway, lactose synthesis and glycerol-3P utilization. For 
lipid metabolism, three new intermediary metabolites 
were introduced (GlycerolSP) and two fatty acid primers 
(Primer (C2:0)CoA and Primer (C4:0)CoA). They corre- 
spond to the primers of two (from acetate) or four (from 
/3-hydroxybutyrate; BHBA) carbon-units necessary to ini- 
tiate fatty acid synthesis before elongation, since in rumi- 
nant mammary glands fatty acids are synthesized from 
acetate or BHBA and not glucose. Other additional reac- 
tions were included to describe milk synthesis, such as 
NADPH synthesis through the isocitrate dehydrogenase 
pathway. The input (net uptake) of long-chain fatty acids 
was not considered in the system since these are not syn- 
thesized by the mammary gland and we hypothesized 
that they are not oxidized within the mammary gland 
of lactating animals, as measured in well-fed goat and 
nonruminant through isotopes [35,37]. This hypothesis 
was sustained by carbon balance measurements based on 
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experimental measures of C02 [28]. This latter hypoth- 
esis was relaxed in a second step (see Figure 3) to study 
how our results are impacted by the introduction of oxi- 
dation of long-chain fatty acid within TCA cycle, as mod- 
eled in alternative contexts [18,20-22]. In addition, the 
action of long-chain fatty acids to regulate enzyme activ- 
ities involved in C(4:0) . . . C(16:0) production was taken 
into account by the fact that C(4:0), . . . C(16:0) fluxes of 
synthesis are explicitly given in the dataset. 

Protein synthesis was summarized by the number of 
peptide links synthesized (peptide synthesis) and the ATP 
required for each peptide synthesis was fixed at 5 in the 
present model [23]. Similarly, protein degradation was 
summarized by the number of hydrolyzed peptide links 
(peptide hydrolysis). The ATP consumption for each pep- 
tide hydrolysis was set at a single ATP [54]. 

Finally, four output reactions were included in the 
model in order to calculate the overall allocations in car- 
bon (R137), nitrogen (7? 138 ), C02 (R139) and Serine (R135). 

As a last step, this model was restricted to the case 
of ruminants. Some reactions, such as fatty acid oxida- 
tions {R20 to R23) and pyruvate synthesis through malic 
enzyme R^) and some input and output fluxes present in 
the generic stoichiometric model [23] were not considered 
to occur in ruminant mammary glands (fructose {R25); 
glucose production from G6P neoglucogenesis (Rw)> ace- 
tone synthesis (#69), BHBA synthesis (R70), acetoacetate 
synthesis (R71), and ureogenesis (i?5o). 

Altogether, this mammary gland model contained 140 
reactions, involving eleven intermediary metabolites, four 
cofactors and ATP, C02, 02 and NH3. The stoichiometry 
of the system implies that the balance (production minus 
consumption) of these four last outputs can be uniquely 
deduced from the choice of a plausible flux distribution. 
The full model is shown in Figure 2 and the corresponding 
stoichiometric model is provided in the supplementary 
webpage (SBML format). 

Additional information: datasets, inputs, outputs, 
additional knowledge, component 

Two datasets [28,29] were used, in relation to the response 
of dairy cow mammary gland metabolism to some 
increases in the protein supply: a control diet (Ctrl) and 
increased protein supply through casein infusion into 
the duodenum (CN). A complementary dataset (HB), 
previously used in a mechanistic model of the mam- 
mary metabolism, was used as a reference response [21]. 
Datasets are shown in Table 1. All data were converted 
into mmol/h/half udder. 

Fluxes of nutrients taken up on a net basis by the mam- 
mary gland were considered as inputs whereas fluxes of 
nutrients produced in milk, such as lactose, milk protein, 
milk fatty acid arranged in triglycerides or C02 pro- 
duced and released in blood, were considered as outputs. 



Special attention was paid to amino acid inputs and out- 
puts, milk protein output, fatty acids synthesized within 
the mammary gland and total triglyceride output accord- 
ing to established calculation rules [21,30]. Milk proteins 
were modeled by their number of peptide links (peptide 
output). The amino acid inputs or outputs corresponded 
to the balance between their real net uptake minus their 
contribution to peptide outputs (utilization in milk pro- 
tein). When the balance of an amino acid was positive, it 
corresponded to an input flux. When it was negative, it 
corresponded to an output flux (see Table 1). Two input 
reactions were also assumed to be inactive in all data 
sets since the corresponding net uptake were not mea- 
sured in the datasets (Table 1): propionate input (R.2%) 
and triglyceride hydrolysis {Res) in long-chain fatty acids 
input. 

Eight additional linear constraints were introduced to 
the fluxes to ensure that the model was relevant bio- 
logically [21]. These are detailed in Tables 4 and 5 and 
allow building matrix B appearing in Eq.(l) f . The rate of 
fatty-acid acylCoA (C(n:m)-acylCoA) was fixed to be three 
times higher than the triglyceride output, assuming that 
100% of milk lipids are triglycerides: it included all fatty 
acid acylation, that is, fatty acid synthesized in the mam- 
mary gland and long fatty acid taken up. 50% of fatty acids 
were assumed to be synthesized from acetate, i.e. from 
primer C(2:0)CoA and the other 50% from BHBA {primer 
C(4:0)CoA), with the exception of C(4:0) fatty acid, which 
was supposed to be synthesized only from BHBA. It was 
assumed that 30% of NADPH was generated by the isoc- 
itrate desydrogenase pathway and 70% of NADPH was 
generated in the pentose phosphate pathway [21]. 

Carbon was the component chosen to compute and 
analyze the allocation of nutrient inputs in nutrient out- 
puts. The stoichiometry of the metabolic network satisfies 
the mass-invariant property according to this component. 
This required a very precise description of C02 produc- 
tion and consumption in the model reactions and was 
validated with external measurements of the C02 balance 
provided with the datasets. 

Availability of supporting data 

The datasets supporting the results of this article are 
included within the article (Table 1). The SBML file 
for the mammary gland model and tools enabling the 
computation of all results are provided in a com- 
panion web page: http://nutritionanalyzer.genouest.org/ 
SupplementaryMaterial. 

Endnotes 

a These functions are: vg + v\9, vg — v\9, V14 + V19, 

VIA ~ V\9, Vis + V W , Vis ~ V\9, V 6i + Vi 9 , V 6i - Vi 9 

b For instance, for the (Crtl) treatment, there is a 
distribution such that vg < 591 and V13 < 266 are lower 
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than the bounds introduced in [26,27], while 

"64/"63 = 0-62 < 0.67 and glucose is the precursor for 

90.79% of the lactose carbon 

c All the computations were performed by using Matlab 
release 2011 on a Xeon E5645 multicore computer. 

d These functions are: v 8 + V19, v 8 + V19 — 50 V15, 

Vg - V19 - V14, V 8 - V19 - Vi 5) V14 + V19, V i4 + V19 - 50 V 15 , 

"14 - "19, "15 + Vi 9) V X5 - Vi 9) V 6 4 + "19, 

"64 + V19 - 50 Vi 5 , V 64 - V i9 - V14, V 64 - V 19 - V 15 . 

e http://nutritionanalyzer.genouest.org/ 
SupplementaryMaterial. The access to the pre-formatted 
model and datasets supporting the results is planned to 
be on request. A valid access is provided by: login : 
SuppMat; password : tamppuss 

'Mathematically, these constraints can be derived as: 
0.7v 82 - 0.3v 58 = 0; v 51 - v 52 = 0; v 54 - v 90 = 0; 
"55 - "91 = 0; v 86 - v 92 = 0; v 87 - v 93 = 0; v 88 - "94 = 0; 
"56 - 3v 62 = 0. 

Additional file 



Additional file 1: SBML version of the metabolic model. The complete 
metabolic model of mammary gland metabolism is provided in the free 
and open standard SBML representation format (Systems Biology Markup 
Language). 
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